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Colloidal systems are often modelled as fluids of hard particles (possibly with an ad- 
ditional soft attraction, e.g. caused by polymers also contained in the suspension). In 
simulations of such systems, the virial theorem cannot be straightforwardly applied to 
obtain the components of the pressure tensor. In systems confined by walls, it is hence 
also not straightforward to extract the excess energy due to the wall (the "wall tension" ) 
H from the pressure tensor anisotropy. A comparative evaluation of several methods to 

^ circumvent this problem is presented, using as examples fluids of hard spheres and the 

Asakura-Oosawa model of colloid-polymer mixtures with a size ratio g = 0.15 (for which 
OO the effect of the polymers can be integrated out to yield an effective attractive potential 

between the colloids). Factors limiting the accuracy of the various methods are carefully 
T~H discussed, and controlling these factors very good mutual agreement between the various 

^ methods is found. 

" ' Keywords: hard spheres; Asakura-Oosawa model; ensemble-mixing method; Monte Carlo 

' ' simulation 



PACS numbers: 82.70.Dd, 05.10.Ln, 61.20.Ja, 68.08.-p 

1. Introduction 

In many colloidal dispersions, the interaction between the colloidal particles (or at 
least the repulsive part of these interactions) can be represented as that of hard par- 
ticles (hard spheres, hard rods, etc.)lil21Ml^ Colloidal dispersions are model systems 
for the experimental study of cooperative phenomena in condensed matter, since 
the size of these particles (in the fj,m range) and their slow dynamics allows their 
study in unprecedented detail, which would not be possible for condensed matter 
systems formed from small molecules. Moreover, hard particles are the archetypical 
many-body systems, for which computer simulation methods were first introduced 
a long time ago: the importance sampling Monte Carlo method ^ was introduced 
discussing a fluid of hard disks; similarly, Molecular Dynamic methods were intro- 
duced in the context of the study of the hard sphere fluid ^1. These fluids formed from 



1 



October 10, 2011 0:16 WSPC/INSTRUCTION FILE 

Methods ' t o ' compute ' Pressure 



2 D. Deb et al. 

hard particles also are an important testing ground for analytical theories of fluids 
The above remarks just sketch a few motivations to explain why the simulation 
of fluids formed from hard particles is of fundamental importance. 

However, simulations of fluids containing hard particles suffer from one impor- 
tant technical difficulty: the well-known virial equation for the pressure tensor (for 
systems where particles interact with pair- wise forces fifij) with Vij — \ fij\,rij = 

p-' = pkBT6'^^ + ±{j:nn,)4)^^^ (1) 

is not straightforwardly applicable (for hard spheres of diameters (7^, the force 
/ = if > (a-i + crj)/2, while distances < ((7^ + o'j)/2 are excluded, and for 
fij = + the potential energy jumps from zero to infinity, and hence the 

force is also not defined in this case). Note that in Eq. ([T]) q,/3 denote Cartesian 
components, p = N/V is the density (7V=particle number, V— "volume" of the 
rf-dimensional system), fc^ is Boltzmann's constant, T absolute temperature, and 
(• • • )nvt is understood as a statistical average in the canonical (NVT) ensemble. 
Computation of the pressure is not only of the interest to clarify the equation of 
state of the considered system in the bulk, but is also required when one considers 
a system confined by walls, and wishes to estimate the excess free energy due to the 
walls. Assuming a geometry where the system is confined by two identical walls at a 
(large) distance D apart, the wall- fluid surface excess free energy per unit area 7^,/ 
(also called "wall tension") is expressed in terms of the anisotropy of the pressure 
tensor as ^ 



In^f = D{Pn - Pt)/2 (2) 

where Pjv — P^^ is the "normal pressure" which in equilibrium is homogeneous 
in the system, while Pt is the average tangential pressure 

D 

Pt = {II2D) J ip-^z)+pyy{z))dz. (3) 


Note that Eq. ^ supposes that the two walls are far enough from each other 
that in between the fluid can attain its bulk properties, i.e. the two walls are non- 
interacting, the density profile p(z) as well as the transverse pressure tensor com- 
ponents p^^(z) and pyy{z) are constant, independent of z, for a broad region of z 
near z = D/2. In contrast, p{z) does depend on z near the walls, as well as p^^{z) 
and pyy{z), while the transverse component p^^{z) is strictly independent of z, irre- 
spective of the choice of D. The generalization of p"^ in Eq.[l] which addresses only 
the behavior in the bulk, to such inhomogeneous systems with external boundaries 
will be discussed in Sec. 3 below. 
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In the present paper, we shaU discuss methods from which the pressure and 
the waU tension can be computed for fluids where particles experience hard-core 
interactions, and hence the straightforward use of Eq. [l] or its extension to inho- 
mogeneous systems is not possible. In Sec. 2, we shall describe the application of a 
recent method due to Schilling and Schmidl^^, to compute absolute free energies of 
disordered structures by molecular simulations, considering the situation described 
above, where a system is confined by walls. We assume planar, structureless walls, 
with a smooth repulsion described by a Weeks-Chandler- Andersen (WCA) l^type 
potential, 

VWcAW=4£[KA)i2-(a„,/z)6 + l/4],0<z<z, = a^2i/6 , (4) 

while Vwca{z > Zc) = 0. Here the strength (e) and range (cr^) of the potential 
provide parameters which can be varied, and so also the wall tension j^f can be 
changed. For the chosen geometry, the free energy density can be written as 

f{p,T,D) = Mp,T) + {2/D)j^f (5) 

where fh{p,T) is the free energy density of a corresponding bulk system at 
density p and temperature T. We shall test the applicability of Eq. ^ for the 
simplest case, namely the hard sphere fluid. 

In Sec. 3 we follow the alternative approach of De Miguel and Jackson 1^21^ where 
Eqs. ([2]), ([3]) are used for hard particles in spite of the fact that a straightforward use 
of Eq. ([1]) is not possible. In this method, virtual changes of simulation box linear 
dimensions are performed to sample the probability that no molecular pair overlaps 
would occur as a result of these virtual moves. From this probability, the pressure 
can be estimated. While we already have verified the accuracy of this method for 
the case of a hard sphere fluid we here extend the method to the case where 
the interparticle potential has both a hard core repulsion and a soft attraction. A 
well-known and practically relevant example for such a situation is provided by the 
Asakura-Oosawa model E l ^^ l ^^ l for colloid-polymer mixtures, in the limit where the 
polymers are much smaller than the colloids. Describing the colloids as hard spheres 
of diameter a and the polymers as soft spheres of diameter (jp, polymer-colloid 
overlap is also strictly forbidden, but polymer coils can interpenetrate and hence 
overlap with negligible energy cost. For q = ap/a < 0.154, one can integrate out 
the polymer degrees of freedom from the partition function exactly, replacing their 
physical effects by an effective attractive interaction between the colloids ISIIEIEZEll 

while VAoif < cr) = 00, Vao{''' > cr -|- CTp) = 0. The strength of this attractive 
potential is controlled by the "polymer reservoir packing fraction" 77^, which is 
related to the chemical potential /ip of the polymers as rf^ = (tt /Q)a'^ exp(pp/kBT). 
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Both methods of Sees 2. and 3 are checked by applying an independent third 
method, extending!^ a thermodynamic integration method where one gradually 
inserts a wall potential in a system without walls (with periodic boundary condi- 
tions). Our extension " called "ensemble mixing method" samples the free energy 
of systems with intermediate Hamiltonians. This new method is presented in Sec. 4, 
while Sec. 5 presents our conclusions. 

2. Estimation of wall tension from computation of absolute free 
energies 

As is well known, the absolute free energy of a system is normally not a straightfor- 
ward output of a Monte Carlo simulation The most commonly used strategy is 
to use thermodynamic integration, supposing that the Hamiltonian 'H(A) depends 
on a parameter A that can be varied from a reference state (characterized by Aq) 
whose free energy is known, to the state of interest (Ai), without crossing any phase 
transition, 

Al 

AF^F(Xi)-F{Xo)^ J dX'{H{X')/dy)x' ■ (7) 

Ao 

For a dense fluid. Schilling and Schmid ^ propose to take as a reference state 
a configuration that is representative for the structure of interest (obtained within 
a standard simulation of the considered system). From this well-equilibrated state 
a reference state is constructed by fixing this particular configuration with suitable 
external potentials that hold all particles rigidly in their position {rf} in that 
particular configuration (and the internal interactions between the particles can be 
switched off). In practice. Schilling and Schmid ^ used as a pinning potential in 
the reference state 

C/,.ef(A) = Xj2^{\n - rP|/re„t), with0(x) =x-l, (8) 

i 

where it is understood that particle i can only be pinned by well i at r/, and not 
by the other wells (but one carries out identity swaps to make the particles indistin- 
guishable). For this non-interacting reference state the (configurational) reference 
free energy is 

Fr,i{X)=HN/V)-Hl + {Vo/V)g4,{X)] , (9) 

3 

where Vq = 47rrj?^j/3 and ^^(A) = 3A^(e'*' — ^ e'^/kl), for the above choice of 

k=0 

(j){x). Then, intermediate models are defined as 
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H'(A) =Hi„t + C/,.ct(A) , (10) 

where Hint describes the interactions in the system which then are switched on 
(if necessary, in several steps). The free energy difference relative to Frcf(A) then is 
computed by thermodynamic integration, for which {d'H,:ci{X)/dX) = (X^i'^d^i ~ 
Tf\/rcut))\ needs to be sampled. For details, how this is done efficiently, we refer to 
the original publication Schilling and Schmid^i^have already tested this method 
for hard spheres in the bulk, obtaining the free energy for several densities, and 
verifying the (expected) agreement with the Carnahan-Starling equation of state 

i2n 

In the present section, we present a straightforward extension of this method 
to a hard sphere system in a. L x L x D geometry, choosing L — 6 (lengths being 
measured in units of a), choosing a packing fraction rj — (7r/6)p = 0.3686 and 
several choices for D (Fig. [ij. In order to obtain a reliable estimation of the error 
of this method, the approach was repeated for each choice of D for 5 independent 
(but equivalent) reference configurations. The total number of Monte Carlo steps 
per particle for each of these runs was in the order of about 4 • 10^ for each of 
the 500 steps of the thermodynamic integration and about 2 • 10^ for each of the 
40 steps in which the potential wells were switched on. Fig. [T] shows that indeed 
F/D plotted vs. 1/D is compatible with a straight line, and the intercept of this 
straight line agrees with the Carnahan-Starling result (which was verified directly 
by also running a L x L x L system at the chosen packing fraction). The accuracy 
with which the straight line can be fitted to the data points allows us to estimate 
from Eq. ^ with a relative accuracy of about 1%, and the result (choosing 
units such that fcsT = 1) 7 = 1.01 ± 0.01 is fully compatible with the result of the 
ensemble mixing method for this case (Sec. 5). Thus, Fig.[l] is a valuable feasibility 
test that shows that the Schilling-Schmid method is useful for the estimation of 
interfacial excess free energies of fluids. 



3. Computation of the pressure tensor in systems of hard particles 
confined by walls 

In this section we consider the AO-model of colloid-polymer mixtures, where par- 
ticles interact both with a hard core interaction (for r < a) and a soft attractive 
part of the interaction, as described in Eq. ([6]). We wish to apply Eq. ([2]), and need 
to compute both P/y £^nd Pt for this purpose. Furthermore we allow again for walls 
where the WCA potential, Eq. Q, acts. 

First we note that in a system which is inhomogeneous in the z-direction the 
generalization of Eq. ([l]) to a local pressure tensor can be written as '221 
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p"^iz) = p{z)kBTS'^' - ^^Y.'-^,{§-)^0[{z zi)lz.mz, - z)/z,,]/\z,,\ 



-p{z)z^^Vwca{z)S'^'S^' . (11) 

Note that in the total force acting on a particle there is also a contribution due 
to the wall potential, which has been written down explicitly in the last term, which 
contributes only for a = z and /3 = z. Thus, we can decompose p^^ into the three 
contributions p{z)kBT and Pwan(^) = —p{z)z{dVwcA{z)/dz) and the part due to 
interparticle interaction p?^^{z), 



p^'^ = p{z)kBT + p^,{z)+p^,u{z) , (12) 
noting that while all three components of p^^ do have a pronounced z- 



dependence in Eq. (12), the total normal pressure p^^ in fact must be independent 
of z. 



In our case Eqs. (Ill, (12 1 can be applied directly with respect to the soft 
attractive part of the interactions, but not with respect to the hard core part. Here 
we apply the principle that all contributions to the pressure simply are additive, 

Pfnt (^) = P"nt, soft (^) + Pfnt, hard (^) ' (13) 

and we only compute soft(-^) directly from the virial expression, i.e. the 
second term on the right hand side of Eq. (jllj), while for Pj"^ hard(-^) apply the 
method of De Miguel and Jackson as already used in our previous work on the 
hard sphere fluid 

In the latter method, one considers virtual volume changes according to 



D^D' = D{1-S),L^L' = L{l-S) (14) 

where ^ << 1. In any equilibrium configuration, there cannot be any overlaps 
between the hard cores of any pair of particles, but in dense fluids these virtual 
moves will create some overlaps. One then can sample the probabability Pnov(0 = 
exp(— &^) w 1 — &^ that no overlap of particles occurs. Denoting the constant h for 
the moves from D to D' as bj^, and for moves from L to L' as bx, we obtain 



Pjv. hard/(pA:sT) = 1 + bN/N, (15) 



where 



-'T, hard 

/{pkBT) = l + bT/N 



(16) 
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D 

PiV.hard = (l/^) /pfnlhardW'^^ : (17) 



and 



D 

PT,hard = (l/2i^)y'[pr„lhardW+Pfn't,hardW]'^^- ■ (18) 


In view of the substantial difficulty in obtaining very accurate simulation "data" , 
we have not attempted to compute the individual profiles p"^^{z), noting that for 
the desired use of Eq. ((2| only the quantities averaged across the film matter, as 



written in Eqs. (17 1, il8\ 



As an example. Fig. ^ shows the density profile and the profile Pwaii(-z) as 

3?ft(z) and pf^^{z) = pyg{ 



well as the profiles of those parts of both p?fj(z) and pfnJ(z) = pfnt(-2) iha.t are due 



to attractive parts of the potential of the AO model (Eq. |6]). We recognize that all 
these contributions exhibit rapid large oscillations near the walls. These oscillations 
need to be well resolved (with a binning of the z-coordinate that is about an order 
of magnitude smaller than the hard sphere diameter, which sets the scale for these 
oscillations) and carefully sampled. Otherwise systematic errors are made when 
these contributions are integrated over z, and despite of the simplification that only 



integrated pressures are needed (cf. Eqs. (17), (18)) all results would suffer from 
systematic errors. 

However, even with this simplification it is very difficult to obtain a very good 
accuracy, because systematic errors that occur if the range of ^ for which PnoviO 
is sampled are not extremely small. 

For the sampling of -Pnov(Oi ^ careful sampling of the probability that pairs of 
particles occur at a distance a + Ar from each other is required, for very small 
values of Ar is required. Of course, the number of particles n{Ar) that occur at a 
distance in between a and a + Ar goes to zero as Ar — > 0. The constants and 



bx in Eqs. (15), (16), as well as the analogous constant for the bulk pressure in a 



system with periodic boundary conditions, depend very sensitively on the magnitude 



of ^ that is considered in the virtual volume changes considered in Eq. (14), or, 
equivalently, the range of Ar that is used to fit the derivative of dn(Ar) / d{Ar) as 
Ar — >■ 0. This problem is illustrated in Fig. |3]for several cases of rjb- Note that the 
range of ^ from < ^ < 5 • 10^"' must be probed to allow an accurate linear fit from 
which the limiting value for ^ — > can be extracted. If one fits a larger range of ^, 
systematic errors do arise, particularly for larger values of the packing fraction. Note 
also that the relative range of pressure variation with ^ increases with increasing 
packing fraction. Thus, it is rather easy to find a correct result for rather small 
values of the bulk packing fraction rji,, but much more difficult when rji, approaches 
the region where the onset of crystallization occurs. This problem is not specific for 
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the AO model, but occurs for the pure hard sphere system as weU. In fact, already 
an examination of the local density near the walls (Fig. [4| shows that in the fluid 
phase both models are still very similar (only the solid phases at the liquid to solid 
transition have rather different packing fractions and the interfacial stiffness 
between the coexisting phases also differ appreciably!^. In the AO model, the first 
peak of the layering structure near the WCA wall is slightly higher and slightly 
sharper as in the pure hard sphere case, but otherwise the differences are rather 
minor. 

Fig. |5] then presents the various contributions to the pressure as well as the wall 
tension 7^,/ as a function of bulk packing fraction. Of course, for the application of 
Eq. ([2]) the individual terms P/v, Pt need to be obtained with an accuracy much 
better than one part in a thousand, if for large D a meaningful accuracy of 7^,/ is 
desired: From Fig. [sja), we recognize that near the crystallization transition (which 
occurs forl^ rjp — 0.494) the pressure in the bulk is of order 10^, while 7^,/ is of 
order 10°. With D/2 of order 20, Eq. ^ implies P/v - Pt is of order 1/20, and 
hence {P]s[ — Pt) / Pn is of order 1/200. This consideration already explains why (on 
the scale of Fig. [sja,) the difference between Pjy and Pt is invisible, and getting 7^,^ 
accurate to a few percent hence is a computational tour de force. This consideration 
explains readily why already for the case of hard spheres it has happened that 
different methods for the estimation of 7^,/ have yielded results that slightly different 
from each other (see Ref . 1^ for a recent account on the literature of this problem) . 

4. The ensemble mixing method 

In Ref. we have already presented very briefly the ensemble mixing method where 
one gradually transforms from a system without walls in the (NVT) ensemble (with 
periodic boundary conditions in all three space directions) to a system with two 
parallel walls at distance D from each other (and periodic boundary conditions 
only in the x^y directions parallel to the walls). Describing the system without 
walls by a Hamiltonian Hi(X), X denoting a point in configuration space, and the 
system with walls by 7i2{X), the idea (originally proposed by Heni and Lowenli^ 
is to construct a smooth path along which one can carry out a thermodynamic 
integration, that yields the free energy difference between the two systems, from 
which hence 7^/ can be inferred. Heni and Lowen applied this idea to the fluid 
of hard spheres, and also the recent extension due to Deb et al. ^ deals with the 
hard sphere fluid only. In the present section, we demonstrate that this method can 
equally well be applied to the AO model. 

Since both systems T-li{X) and H2 {X) refer to the same choice of particle number 
N and volume V, i.e. the same phase space, one can think of constructing an 
interpolation Hamiltonian T-L{X) as follows, 

H{X) = {\ - k)Hi{X) + >,:H2{,X) (19) 
where the parameter n G [0, 1] is varied for calculating the free energy difference 
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between the systems (1,2). In our simulations, k is discretized in a set of discrete 
values {Ki}, and the system is allowed to jump from Hi to a neighboring value, 
Ki+i or These moves are accepted or rejected with a Metropolis step, and 

hence in order to ensure a large enough acceptance probability, the number of 
intermediate systems {k^} needs to be chosen proportional to the volume V of 
the system. We sample the relative probability P{i) to reside in state i with a 
variant of Wang-Landau samplingUHlH Alternatively, other free-energy schemes like 
successive umbrella sampling can be appliecl^. The free energy difference between 
the two states i and i -I- 1 is given by kBT[lnP{i) — \nP{i + 1)]. In this way one 
can sample the free energy difference AF(D) for any given choice of D, between a 
system with walls and the system of corresponding thickness and particle number 
but without walls. The wall free energy then follows as 



7^/ = lim AF{D)/{2L^kBT) . (20) 

Note that due to the excess density ps at the walls (which is present in the 
systems with walls, so that p = N/V — pb + {2/D)ps), for any finite D the bulk 
density of the system with periodic boundary conditions (which is p) differs slightly 
but systematically from the bulk density in the system with the walls (which is 
Pb = p — {2/D)ps). This fact already suggests that the extrapolation in Eq. (20) 
should be carried out in terms of a plot of AF{D)/{2L'^kBT) linearly versus 1/D 
(Fig. |6]). In this way the data included already in Fig. [5](b) have been generated. 

Of course, in this method the lateral linear dimension (L = S.OScr) is rather small, 
and often even smaller choices for L have been used (e.g., L = 5a)^. However, no 
significant finite size effects associated with the use of such small values of L were 
detected. For the pressure tensor calculations, somewhat larger choices for L have 
been possible (near L = 12). Again, the good agreement between both methods 
to obtain j^,f for a significant range of rjb (Fig. ^jp) can be taken as an indication 
that neither method suffers from significant systematic errors, in the range of the 
parameters of interest for the present model. 



5. Conclusion 

In this paper, a comparative study of several methods to compute the excess free 
energy of walls in fiuids (the "fiuid-wall tension") has been presented, considering 
fluids where the particles exhibit a hard-core repulsion, and hence the standard 
method based on the wall-induced anisotropy of the pressure tensor is not straight- 
forward to apply. We have described three methods: (i) In the Schmid-Schilling 
method, the absolute free energy of the system is computed. When one does this 
for systems in a film geometry with several choices for the distance D between 
the walls confining the film, one can infer the wall tension from the coefficient 
of the term in the free energy that exhibits a 1/D variation. We have tested this 
method so far for the hard sphere fluid only, but we did find good agreement with 
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other methods, (ii) The pressure tensor anisotropy method can in fact be gener- 
ahzed to systems, with hard potentials, by a careful sampling of the probability 
Pnov(C) that particles do not overlap when virtual changes of the linear dimensions 
by a factor 1 — ^ are performed. We have applied this method to the AO model of 
a colloid-polymer mixture, where the (implicitly treated) polymers are responsible 
for an effective attraction between the particles. The latter part of the potential 
is treated via the virial formula for the pressure tensor, while the repulsive part is 
treated with a sampling of Pnov(C)i mentioned above. We have shown that by a 
careful consideration of accuracy issues this method does yield valid results, but a 
substantial computational effort is needed, (iii) The third method is the ensemble 
mixing method, where one considers systems with periodic boundary conditions and 
no walls "mixed" (in the sense of a combined Hamiltonian, Eq. ( 19 )) with a system 
with the same particle number and linear dimensions, but bounded by walls. This 
method has been used to test the other two methods, and our practical experience 
is that this last method amounts to the relatively smallest computational effort, in 
comparison with the other two methods. In the future work, we plan to use such 
methods to study the wetting behavior of colloidal dispersions. 
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Fig. 1. Extrapolation of the free energy /(»?c, D) of a hard sphere system at a packing fraction 
ric = 0.3686 of colloids plotted versus 2/D, D being the distance between two parallel walls, 
at which the WCA potential (Eq. |4|) acts. Note that the intercept at the ordinate is known 
accurately from the Carnahan-Starling equation of state. The slope of the straight line fit yields 
the wall tension. 




Fig. 2. Profiles of the packing fraction r](z) and wall contribution P^^f^{z) to the normal pressure 
Pwall(z), case (a), and of the attractive potential parts to both normal {Pg^fti^)} and tangential 
'tPsoft(^) = PgQft(2)} pressure tensor, for the AO model at bulk packing fraction rji, = 0.4785133. 
Note that pf^f^Cz), Paoft(^) are obtained from the second term on the right hand side of Eq. ([ij. 
Insets show the variations near the left wall on magnified abscissa scales. 
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Fig. 3. Extrapolation of the repulsive part of the total normal pressure Pjv hard a-nd the total 
tangential pressure Pt, hard resulting from the hard sphere interactions at three packing fractions 
plotted vs g for % = 0.2786226(a), 0.3539615(b) and 0.4538897(c). 
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Fig. 4. Density profile p{z) versus z near WCA walls with e = 1 at z = 0, for the hard sphere 
model (a) and the AO model (b). 
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Fig. 5. (a) Various contributions to the normal pressure Pjy (i'jv.walli ^iv.int) and the tangential 
pressure (Pt) for the AO model plotted versus rji,, choosing D = 44.55737 and two WCA walls 
with e = 1.0. (b) Wall tension 7,„j plotted vs. rji,. Open squares are due the ensemble mixing 
method, full squares result from the application of Eq. 
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Fig. 6. Typical example showing how wall tensions 'y^f are extracted using the ensemble mixing 
method for the AO model. Part (a) shows the free energy difference AF{k) / kgT vs. k for e = 2, 
r)i, = 0.4697. Part (b) shows the extrapolation to 1/D —> for several rj^^ 
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